Vertex-centroid finite volume sclieme on tetraliedral grids for 

conservation laws 

Praveen Chandrashekar, Ashish Garg 
TIFR Center for Applicable Mathematics, Bangalore-560065, India 



Abstract 

Vertex-centroid schemes are cell-centered finite volume schemes for conservation laws 
which make use of vertex values to construct high resolution schemes. The vertex values 
must be obtained through a consistent averaging (interpolation) procedure. A modified 
interpolation scheme is proposed which is better than existing schemes in giving positive 
weights in the interpolation formula. A simplified reconstruction scheme is also proposed 
which is also more accurate and efficient. For scalar conservation laws, we develop limited 
versions of the schemes which are stable in maximum norm by constructing suitable 
limiters. The schemes are applied to compressible flows governed by the Euler equations 
of inviscid gas dynamics. 

Keywords: Finite volume method, unstructured grids, reconstruction, maximum 
principle, compressible flows 



1. Introduction 

Finite volume methods on unstructured grids have become very useful in solving con- 
servation laws, especially those arising in fluid dynamics and in which the computational 
domain is of a complex shape. This is mainly due to the maturity of grid generation 
tools, especially triangular and tetrahedral grids. Also, the local conservation property 
of finite volume methods allows them to automatically compute discontinuous solutions. 
Finite volume methods can be classified based on the location where the solution values 
are stored [21 [H]. In cell-centered methods, the solution is assumed to be stored at the 
centroid of the cells, while in vertex- centered methods, the solution is stored at the vertices 
of the grid. In the latter case, a cell has to be constructed around each vertex for which 
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there are several possible alternatives. The most commonly used method is the dual grid 
scheme in which the cell is obtained by joining the cell centroids to the face and edge 
centroids. Strictly speaking, the unknowns in a finite volume method are cell averages, 
which are not associated with any particular spatial location in the cell. However, the 
cell average value gives a second order accurate approximation to the solution at the cell 
centroid. For constructing second order accurate schemes, the cell average value can be 
taken to be the point value at the cell center, which is the approach taken in the present 
work. 

The construction of second order accurate schemes leads to several possible alterna- 
tives. Since our work uses cell centered schemes, we restrict the discussion to this case. 
The finite volume method updates the cell average value of the solution, and the detailed 
variation of the solution inside the cell is not known. To construct higher order schemes, 
the solution inside each cell must be reconstructed by using the cell average values of the 
neigbouring cells. The common approach is to reconstruct the solution variation in each 
cell by a linear polynomial by making use of the cell average values in the current cell 
and some of its neighbours, which forms the stencil of the reconstruction. This can be 
achieved by first approximating the derivatives of the solution in the cell by using either 
a Green-Gauss procedure or a least squares procedure [2]. The reconstructed solution 
may be limited in order to maintain oscillation-free solutions. While in one dimensional 
problems, the TVD criterion [8] can be used, for multi-dimensional problems on unstruc- 
tured grids, the usual criterion is to ensure that the reconstructed solution lies within the 
minimum and maximum of the neigbouring cell average values. The two states obtained 
on either side of a cell face by the reconstructed solution are used to compute the fiux 
across the face using some numerical fiux function. 

The vertex-centroid scheme [U |5] is a finite volume scheme for conservation laws on 
triangles and tetrahedral grids. The basic unknowns are the cell center values but it also 
makes use of vertex values in order to construct the second order scheme. The authors 
in [U [S] obtain a formula for the reconstructed value at the center of a cell face by making 
use of the geometrical properties of triangles/tetrahedra, and using the cell center and 
vertex values of the solution. There is no need to explicitly compute the derivatives of 
the solution in each cell as in the usual cell-centered finite volume schemes. The vertex 
values are obtained through an interpolation/averaging procedure that is exact for linear 
polynomials. This procedure has been refered to as pseudo-laplacian averaging in the 
literature. The interpolation formula for each vertex is a linear combination of the cell- 
center values adjacent to that vertex and it is desirable from a stability point of view that 
the weights in this formula should be positive. Applications to compressible transonic 
flows with weak shocks have been given in |1[ E] and the vertex-centroid scheme was able 
to compute stable and accurate solutions even without a limiter, demonstrating good 
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stability properties of the scheme. For stronger shocks, the authors recommend using a 
minmod hmiter though a stabihty analysis was not shown. 

A vertex-centroid scheme which is similar to the schemes in [U |5] has been used in 
for two dimensional compressible flows, but it uses an area averaging procedure which 
is not second order accurate on general meshes. However the authors use very regular 
meshes and find that the solutions have extremely low levels of numerical dissipation. 
They also comment that the vertex-centroid scheme simplified the inter-mesh transfers 
of an unstructured mesh multi-grid scheme, in both the fine-to-coarse and coarse-to-fine 
directions. 

In this paper, we analyze and improve the vertex-centroid schemes in three respects. 
We present a modified averaging procedure which reduces the number of negative weights. 
With the new scheme the matrices arising in the determination of the weights have better 
behaviour. Secondly, we analyze the accuracy of Frink reconstruction procedure and also 
propose a simplied procedure which is theoretically more accurate. Thirdly, we develop 
limited versions of the schemes which can be shown to be stable in maximum norm when 
applied to scalar conservation laws. Under a CFL condition, the solution at the new time 
level in each cell is shown to be bounded between the solution at the neigbouring cells 
and the average values at the vertices of the cell. This is an important design criterion for 
numerical schemes solving hyperbolic conservation laws since they can develop discontin- 
uous solutions, even starting from smooth initial data. The developed schemes are tested 
on several problems governed by the inviscid Euler equations modeling compressible flows 
including those with discontinuous solutions. The numerical solutions are compared with 
exact solutions or with experimental data to demonstrate the performance of the schemes. 

The rest of the paper is organized as follows. The basic formulation of the vertex- 
centroid scheme is explained. The vertex interpolation schemes are discussed and a new 
interpolation scheme is proposed. We then study the error in the reconstruction schemes 
of Frink and a simpler reconstruction scheme proposed here which we refer to as upwind 
reconstruction. Limited versions of both the reconstruction schemes are described and 
their maximum stability is analyzed for a scalar conservation law in two and three dimen- 
sions. Finally, the schemes are applied to problems governed by Euler equations of gas 
dynamics and their performance in computing discontinuous solutions is demonstrated. 

2. Vertex-centroid scheme 

Consider a system of conservation laws in d space dimensions which can be written as 




(1) 
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Figure 1: Reconstruction in vertex-centroid scheme 



Here U is the vector of conserved variables and Fj, i = 1, . . . , are the Cartesian compo- 
nents of the flux vector and d is the number of spatial dimensions. We solve this equation 
numerically on a triangulation made of triangles in 2-D and tetrahedra in 3-D. We assume 
that the cells and vertices are numbered in some way; the z'th cell is denoted by Cj. The 
face between the i'th cell and the j'th cell is denoted Sij. The set of indices of the cells 
sharing a face with Ci is denoted by N{i). Then the semi-discrete finite volume scheme 
for cell Ci is 

E H{Ut,.Ur-^.n,^) = ^ (2) 

where riij G is the normal to Sij and pointing into cell Cj with \nij\ = \Sij\, | ■ | denotes 
the measure of the corresponding geometric objects and if is a numerical flux function. 
Note that we have used a mid-point quadrature rule to approximate the flux integral 
on the cell faces which is sufficient to obtain second order accuracy. The numerical flux 
function must satisfy the following conditions: 

• Consistency 

d 

H{U,U,n) = Y,Fi{U)ni 

i=l 

• Conservation 

H{Ui,U2,n) + H{U2,Ui,~n) = 

In the scalar conservation law case, the numerical flux function is a monotone; in this 
case, the numerical flux H{a, b, n) is an increasing function of the first argument and a 
decreasing function of the second argument. 
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The values U^^ and are the two states at the center of face Sij obtained by some 
reconstruction process in the cells Ci and Cj respectively, see figure m. In any tri- 
angle/tetrahedron, the line joining any vertex to the center of the opposite face pasess 
through the centroid of the cell. The ratio of the distance between any vertex and the cell 
centroid, to the distance between the cell centroid and the center of the opposite face is 
2:1 for a triangle and 3:1 for a tetrahedron; these properties are used in the reconstruction 
schemes. Let Vij be the value at the vertex of cell Ci which is opposite to the face Sij 
and Wij denote the arithmetic average of the vertex values on the face Sij. Note that 
Wij = Wji since both the values correspond to the same face, but Vij ^ Vji since they 
refer to different vertices. Then the reconstructed values on either side of a face Sij is 
given by Frink |3j as 

U+ = U, + aaiWij - V,j), Ur^ = U, + a,{W,j - V,i) (3) 

where the subscript d refers to the spatial dimension; for 2-D we have a2 = 1/3 while in 
3-D, we have 0:3 = 1/4. We will refer to this as the Frink scheme. 

An alternate reconstruction which does not make use of face averages is given by 

U^ = Ui + /3d{Ui - Vj), Ur] = U, + PdiU, - V,i) (4) 

where [32 = 1/2 and /^s = 1/3. We will refer to this as the upwind scheme. This scheme 
is a simple consequence of extrapolating the solution at the cell center and vertex to the 
center of the opposite face. Since this scheme uses states on only one side of the face for the 
reconstruction, we refer to this as the upwind reconstruction scheme. Computationally, 
the scheme given by equation ^ is more simple and efficient compared to ^ since it 
does not require the face average values Wij. Later we will see that the upwind scheme 
is also advantageous when limiters are used, since it leads to less restrictions on the 
reconstruction to be stable in maximum norm. 

3. Vertex interpolation scheme 

The basic finite volume scheme used here is based on cell center values which are 
updated by the finite volume scheme. In this case, we do not know the solution at the 
vertices of the grid, which has to be approximated through some consistent interpolation 
procedure. To discuss the interpolation scheme, we take a generic vertex and assume 
without loss of generality, that it is located at the origin. We then find all the cells 
which contain this vertex; let the centroids of these cells have coordinates {xj,yj, Zj), 
j = 1, . . . ,n and let Vj be the Euclidean distance from the j'th centroid to the vertex. 
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In the literature, the following interpolation procedures have been used, an area/volume 
averaging procedure and an inverse distance averaging procedure, given by 



V = — or V = — 

where the summations are over all the n cells containing the vertex. The inverse distance 
averaging can be viewed as a local Shepard interpolation [18]. On general meshes, both 
the above procedures are only first order accurate, i.e., they are exact only for constant 
functions. On uniform and isotropic grids, they become second order accurate. Note that 
for the finite volume scheme to be second order accurate the interpolation scheme must be 
atleast second order accurate. In [10], the area averaging procedure has been used together 
with high quality meshes leading to very accurate solutions. While these methods may 
not be accurate in general, they are very robust since they are convex combinations of 
the cell center values and hence the interpolated value is bounded between the minimum 
and maximum of the cell center values. This is a useful property while solving hyperbolic 
conservation laws which can have discontinuous solutions. 

3.1. Laplacian averaging 

In order to construct an interpolation procedure which is exact for linear polynomi- 
als, Frink [1] extended a method for 2-D Navier-Stokes solutions given in [9j to three 
dimensions. In this approach, we first assume an interpolation formula of the type 

V = (5 

The necessary conditions for this formula to be exact for linear polynomials are 

WjXj = 0, WjVj = 0, Y ^i^i = (6) 
j j j 

In order to determine the weights Wj, the following minimization problem is solved with 
respect to the weights together with the constraints given by equations (|6| 



1 

This approach tries to keep the weights as close to unity as possible while still satisfy- 
ing the consistency conditions. In order to solve the constrained minimization problem. 
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we introduce Lagrange multipliers Aa;,Ay,A^; the constrained minimization problem is 
converted into an unconstrained problem in which we have to now minimize the function 



WjZj 



with respect to the weights {wj} and the Lagrange multipliers X^^Xy^X-^ as independent 
variables. The solution is obtained by setting the derivative of the above function with 
respect to each of these variables to zero; the weights are then given by 



where the Lagrange multipliers are obtained by solving the following system of equations 





Ax 








Xy 















(8) 



In |1] the above set of equations is solved explictly using Cramers rule. 

3.2. Consistent Shepard interpolation 

We propose a modification of the Shepard interpolation procedure which makes it 
exact for linear polynomials. In this case, we begin by assuming an interpolation formula 
of the form 



V 



E 



jr. 3 



E 



where the weights should now satisfy the consistency conditions 



(9) 



(10) 



In order to determine the weights Wj, the minimization probem J7| is solved with respect 
to the weights together with the constraints given by equations (10). In order to solve the 
constrained minimization, we again introduce Lagrange multipliers A^,, Ay, A^ in terms of 
which the solution can be written as 



1 + A^ h A^ h A^ — 
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where the Lagrange multiphers are obtained by solving the following system of equations 



E 



S (?) (^ 

e(s)^ 

S(?) (^ 



E 
E 



E 







E 


\y 




E 


K 




E 



(11) 



We solve the above matrix equation explicitly using Cramer's rule. 

Remarks. A first observation in comparing the above two methods is that all the formulae 
in the modified Shepard interpolation are independent of the spatial scales. If /i is a 
measure of the local grid size, then the determinant of the matrix in equation ^ is 0{h^) 
and the numerical value of this determinant can become very small. For example, if 
h = 10~^, then the determinant could be O(10~^^). In the examples, we find that the 
determinant can reach the levels of machine precision even for inviscid grids. Due to 
this reason, the solution of equation (jsj) will also suffer from round-off errors. On the 



other hand, the determinant of the matrix in equation (11) is independent of h and the 
determinant is found to be well behaved, as given in the later examples. 

Both of the above methods do not guarantee that the weights wj will be positive. A 
necessary condition for equations ^ or (10) to be satisfied with positive weights is that 
the cell centers must be distributed on all sides of the vertex. This means that there must 
be atleast one cell with Xj < and atleast one cell with xj > 0, with similar conditions 
for the other two coordinates. Also, these conditions must be satisfied for any orientation 
of the coordinate axes. In equations ^ the magnitude of the different terms depends on 
the spatial distribution of the cell centers around the vertex while in (10) the magnitudes 



depend only on the angular distribution of the cell centers. We conjecture that due to this 
reason, the modified Shepard interpolation scheme will be able to preserve the positivity 
of the weights better than the Laplacian averaging scheme. We indeed find this to be 
the case in our numerical tests given in the later sections. In many cases, the modified 
scheme gives all positive weights while the Laplacian scheme yields some negative weights. 
In other cases, the number of vertices with negative weights is considerably reduced with 
the modified scheme. 

Due to the presence of the inverse distance factor in equation ([o]), the modified Shepard 
interpolation scheme gives more weightage to the cells which are closer to the given vertex. 
This might be beneficial for anisotropic grids which would arise if we perform anisotropic 
grid adaptation, or in the case of boundary layer grids. However in the present work, we 
have only considered inviscid problems for which the grids are not highly anistropic. 
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4. Accuracy of reconstruction schemes 

The derivation of the reconstruction scheme ([s]) makes use of the geometrical properties 
of triangles and tetrahedra, leading to second order accuracy. This scheme makes use of 
all the vertex values including the cell center value. A simpler scheme which is also second 
order accurate is given by equation Q. This is based on a one dimensional extrapolation 
of the vertex and cell center values to the mid-point of the opposite face. In this section 
we study the accuracy of these two schemes by assuming that the vertex values are exact. 

Assume that U is twice continuously differentiable function of the space coordinates. 
We denote by D^U the symmetric matrix of second derivatives. In any cell C, we have 

\\D^U\\ < K (12) 

for some positive constant where the matrix norm || ■ || is induced by the I2 norm for 
vectors [TIJ. 

4-1 ■ 2-D reconstruction 

Consider a triangle C whose vertex values are Vi, V2, V3 and whose cell center value is 
t/o, see figure (|2]). Without loss of generality, assume that the cell center is at the origin; 
then the position vector of the vertices ri, r2, G M? are such that ri + r2 + = 0. Let 
us look at the reconstructed value on the face formed by vertices 2 and 3. The position 
vector of the center of this face is \{r2 + r^) = —\ri. The exact value at the face center 
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is given by a Taylor series with remainder (TU] as 

Ue = Uo- ^r^Go + ^r'^Ho.n (13) 

where we denote Go = Vt/o and Hqi = D^f/(— sri/2) for some s G [0, 1]. Using another 
Taylor series with remainder term, we can write the vertex values in terms of the cell 
center values as 

' ^" (14) 



Uo + rjGo + -rjH,r, 



where Hj = D^U{sjrj) e for some Sj G [0, 1] and j = 1, 2, 3. 
4-1.1. Frink scheme 

The reconstructed value on the face formed by vertices 2, 3 is according to equation ^ 



U+ = Un 



-(V2 + V3)-Vi 



which upon using equation ( 14 ) leads to 

1 



Using equation (13, the error in the interpolated value can be bounded as 

1 



7 



K\\ri\ 



24 " " 12 
where the diameter of the cell is defined as 



II ill ^2 II ^11 - 24 



(15) 



h = max llrJI 

l<i<3 

4.1.2. Upwind scheme 

The reconstructed value is given by equation Q as 



U'- = Uo + - [Uo - V,] 



which upon using equation (14) leads to 



The error in this approximation can be bounded as 



Ue\<-K\\r,\\^ <-Kh'' 



(16) 



We see that the upwind scheme given by equation (16) has smaller error constant than 
the Frink scheme given by equation (15). 
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4-2. 3-D reconstruction 

Consider a tetrahedron C whose vertex values are Vi, V2, V3, V4 and whose cell center 
value is Uq, see figure Without loss of generality, assume that the cell center is 
at the origin; then the position vectors of the vertices ri,r2,r3,r4 G M.^ are such that 
ri+r2-\-rs-\-r4 = 0. Let us look at the reconstructed value on the face formed by vertices 
2, 3 and 4. The position vector of the center of this face is |(r2 + 7*3 + r^) = —^ri. The 
exact value at the face center is given by a Taylor expansion with remainder term as 

U, = U^- =:riGn + ^riH^.r, (17) 



where we denote Go = Vt/o and Hqi = D'^U{—sri/3) for some s G [0, 1]. Using another 
Taylor series with remainder term, we can write the vertex values in terms of the cell 
center values as 



18) 



where Hj = D'^U^SjTj) G M?^^ for some Sj G [0, 1] and j 



1,...,4. 



4-2.1. Frink scheme 

The reconstructed value on the face formed by vertices 2, 3 and 4 according to equa- 
tion (g IS 



'iV2 + Vs + V,)-V, 



(19) 



which upon using equation (18) leads to 

1 



Using equation (17), the error in the interpolated value can be bounded as 

13 



\u- 



el _ ^2 II 111 ^ 24 " ^" 24 " ^" 24 " ^" - 36 



where the cell diameter is defined as 



h 



max \\rj 
i<i<4 



(20) 
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4-2.2. Upwind scheme 

The reconstructed value is given by equation Q as 



(21) 



which upon using equation (18) leads to 



6 O 



Using equation (17), the error in this approximation can be bounded as 

I el _ g II 111 _ g 



We see that the upwind scheme given by equation (21) has smaller error constant than 



the Frink scheme given by equation (19). 



Remark. This analysis shows that the simpler upwind scheme is actually slightly more 
accurate compared to the Frink scheme. However both methods are second order accurate. 
The difference in the error constant does not translate into any appreciable difference in 
the numerical solutions obtained by the two methods, atleast for the test cases that we 
have studied. The purpose of this analysis was to show that there is no advantage in 
using the more elaborate scheme given by equation ([s]) and the more simpler scheme of 
equation (|4]) is more accurate in theory. In the next section where we construct schemes 
stable in the maximum norm, we find that the upwind scheme is more beneficial since it 
requires less restrictions to satisfy the maximum stability. 



5. Maximum stability for scalar conservation law 

Solutions of nonlinear hyperbolic PDEs can develop discontinuous even if the initial 
condition is smooth. While first order accurate finite volume methods on unstructured 
grids are capable of computing discontinuous solutions in a stable manner [3], the higher 
order numerical schemes might produce spurious oscillations near discontinuities. The 
classic approach is to construct total variation diminishing schemes which will prevent 
the generation of oscillations. The TVD concept however is difficult to extend to multiple 
dimensions and unstructured grids. Moreover, there is a negative result by Goodman and 
LeVeque |7] that a TVD scheme can be at most first order accurate in more than one 
spatial dimension. Due to these reasons, one constructs schemes which are stable in the 
maximum norm. This prevents the generation of new extrema and hence oscillations in 
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the solution will be avoided. The first order accurate finite volume schemes combined 
with a monotone fiux are stable in the maximum norm. The basic idea to achieve this 
stability for higher order schemes is to reduce the accuracy of reconstruction to first order 
wherever the solution has a discontinuity and might lead to oscillations. 

In the case of structured grids, there are well established conditions that can be used 
to check the positivity of the schemes. For one dimensional probems, we can write the 
second order semi-discrete scheme as 

= A{u,+i - a,) + - u,) 

If all the coefficients A, B are positive, then the scheme will be stable in the maximum 
norm under a CFL condition. In particular, if Ui is a local maximum, then it will not 
increase in time, and similarly, if Ui is a local minimum, it will not decrease in time. For 
unstructured grids, the cells are not arranged along grid lines which makes it difficult to 
write the scheme in the above form. However, for the vertex-centroid scheme, this can 
be achieved by using the vertex values along with the cell-center values, which is the key 
step in the present construction of limited schemes. 

We propose to write the semi-discrete vertex-centroid finite volume scheme in the 
following form, which not only contains the difference of cell average values but also the 
difference between vertex value and cell average value. 

|C,|^ = H \MUj - + B,j{V,, - (22) 

We discretize the above set of ordinary differential equations using Euler time stepping 
scheme which leads to the following fully discrete update equation 



ra+l 



^" + ^ E "^^^U, + ^ E B^,V,^ (23) 



If all the coefficients A, B in the above equation are non-negative, then the resulting 
scheme satisfies a maximum principle under a CFL condition, i.e., 

min {f/f , Vi;} < < max {f/f , f/;, Vi;} (24) 

The vertex values are obtained by an interpolation formula of the type given in equation 
(g or (g. If all the weights in the interpolation formula are positive, then the vertex 
values are bounded by the cell average values. The higher order accurate scheme must 
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be constructed to ensure the positivity of the coefficients in equation (22) so that spu- 
rious oscillations near shocks are prevented. If the second order reconstruction given by 
equation ^ or equation Q is written as 

U^ = Ui + , U-^ = U, + Af/,, (25) 

then the corresponding limited reconstruction scheme is defined to be 

U± = Ui + %A?7,„ U-^ = U, + %At/,, (26) 

where the limiter 9 G [0, 1] should be chosen to satisy the positivity conditions. When 
6 = 1, the scheme is second order accurate while if 6' = 0, it becomes ffist order accurate. 
Consider a numerical flux function which can be written as 

H{a,b,n) = H+{a,n) + H-{b,n) (27) 

with being a non-decreasing function and H being a non-increasing function. Then 
from the consistency of the numerical flux function, 

j&N{i) 

we can separate the positive and negative fluxes and write equation ^ as 



jeN{i) 



- 2^ rr+ _ rr ^ijAUij 









P^J>0 


H 




Uij) - H'{Ui,nij) 






V ' 



j&N{i)^ 



Due to the properties of the flux function, we know the sign of the terms defined as Pij 
and Qij. We next derive conditions on the limiter function 6 for the two reconstruction 
schemes so that when they are written in the form of equation (22 ), they have non-negative 
coefficients. 
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5.1. Frink scheme 

For the reconstruction scheme given by equation ([s]), i.e., Af/jj = UdiWij — Vij), the 
semi-discrete scheme can be written as 



E 



{21 



Defining 



adiW.j - Vi 



The above scheme can be written as 



\a 



' dt 



j€N{i) 



Pij^ij^ij-rr TT^V Ui) -\- Qij 



^ 2 



{U, - Ui) (29) 



The first term contains the reconstruction in cell Cj while the second term is due to the 
reconstruction in cell Cj. From the coefficients of the terms containing Pij and Qjj, this 
scheme will have positive coefficients provided the following two conditions are satisfied 
by the limiter function 6'jj, 



< 0, 

' 'V,,,-Ui - ' 



1 



Oifij > 



The above two conditions will be satisfied if we choose the limiter as follows 



{U,-Ui){U,~V,)<Q 
{vij ) otherwise 



(30) 



where the function 6{r) must satisfy the following conditions: 

1. e{r) = for all r < 

2. < e{r) < 1 for all r G [0, oo) 

3. 0{r) < I for all r > 

These conditions are satisfied if we choose 

Q(r) = max(0, min(l, 2/r)) 



(31) 



The function ^(r) is sketched in figure ([3]). The first condition in equation (30) checks 
whether Ui lies between the values Vij and Uj] if it does not lie in this range, the limiter 
is set to zero and the order of the reconstruction is reduced to first order. 
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Figure 3: Limiter function 0{r) 
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5.2. Upwind scheme 

For the reconstruction scheme given by equation Q, i.e., A[/jj = (3d{Ui 



Vij), the 



semi-discrete scheme can be written as 



jeN(i) 



1 + 



(32) 



In this case, the first term inside the curly braces aheady satisfies the positivity condition. 
The coefficient in the second term is positive if we choose 6ij = 6{rij), with 6{r) being 
given by equation (31) and 

'''' \{U,~U,) l{U,-Ui) 
We note in this case that there are less conditions on the limiter function 6ij as compared 



to the previous scheme since we do not have to enforce the first condition in equation (30 ) 



This is due to the fact that the reconstruction scheme does not make use of face average 
values Wij but only the vertex and cell values. 

6. Numerical results for Euler equations 

We apply the schemes developed here to inviscid compressible fiows governed by Eu- 
ler equations of gas dynamics. These are a system of hyperbolic PDE representing the 
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conservation of mass, momentum and energy. In three space dimensions, the vector of 
unknowns U and the Cartesian components of the flux vector are given by 



U 



' p ' 




pUi 


pUi 




p6ii + puiUi 


PU2 




p6i2 + pU2Ui 


PU3 




pSis + pu^Ui 


E 




_ {E + p)ui _ 



(33) 



where p is the density, {ui,U2,U3) are the Cartesian components of the velocity, p is the 
pressure and E is the total energy per unit volume given by 



E 



pe + -p\u\ 



(34) 



which is the sum of internal and kinetic energy. For an ideal gas, the pressure is related 
to the density and internal energy per unit mass e by 



P = (7 - 1)P^ 



(35) 



where 7 is the ratio of specific heats and is a constant for a given gas. In all the numerical 
examples, we take 7 = 1.4 which corresponds to air under normal conditions. 

The Euler equations are solved on tetrahedral grids using a matrix-free LUSGS scheme [T7] . 
In the case of unsteady problems, the 3-stage strong stability preserving Runge-Kutta 
scheme is used [19]. As in [3], we reconstruct the primitive variables {p,ui,U2,U3,p) 
rather than the conserved variables. This is beneficial in maintaining the positivity of 
density and pressure. The numerical flux function is either based on the Roe scheme [15] 
or the KFVS scheme [13]. The latter scheme is known to be entropy consistent [12] and 
hence very robust; we use it for the high Mach number problems which contain strong 
shocks. Most of the grids used in this work were generated using the open source tool 

GMSH [g. 

We also compare the current schemes with the limited schemes of Jameson [TO] where 
a central limiter is used. The reconstructed values on any face Stj are given by 



Ut^=U. + -L{U,-V.„V,.-Uj), 



1 



u^=u,--m-v,,,v,.-u,] 



(36) 



where L(a, b) is a limited average which is taken to be 



L{a,b) = -{a + b)[l- R{a,b)], R{a,b) 



a — b 



max(|a| + e/i^/^) 



ge[l,3] (37) 
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Case 


Pi 


Ul 


Pi 


Pr 




Pr 


Test 1 


1 





1 


0.125 





0.1 


Test 2 


1 


-2 


0.4 


1 


2 


0.4 



Table 1: Initial conditions for shock tube problem 



If g = 1, then the above hmiter reduces to the minmod hmiter. Note that larger values of 
q lead to less limiting. For this scheme, we are not able to prove the maximum stability 
property. Note that the term Jameson scheme is used to refer to the reconstruction 



scheme given by equations (36) and (37). 

All the results are obtained using the modified Shepard interpolation. For some prob- 
lems, a few of the vertices have negative weights which we do not modify. In each case, we 
indicate the number of negative weights using the original interpolation scheme of Frink 
and the modified scheme. 

6.1. Shock tube problem 

The standard shock tube problem is solved on a three dimensional grid. The compu- 
tational domain is in the form of a channel with a square cross-section and whose axis 
is along the x-axis. The side walls of the channel are treated as slip walls. The com- 
putational grid contains about 100 points along the axis of the channel. The solution is 
taken along the center line of the channel for comparison with the exact solutions of the 
Riemann problem. Two standard Riemann problems are considered, the Sod problem 
and a problem which involves low densities in the solution. In both cases, the initial 
discontinuity is placed at x = 0.5. The initial conditions defining the Riemann problem 
are shown in table ([TJ. For Test 1 we use the Roe flux function while for Test 2, we 
use the KFVS flux. Limited versions of the reconstruction schemes are used since the 
solutions have shock and contact discontinuities. With the original interpolation scheme, 
91 vertices have some negative weights with the smallest weight being -0.03, while the 
smallest determinant is 1.44 x 10~^^. With the modifled interpolation scheme, there is 
only one vertex with a negative weight of —1.5 x 10^^ while the smallest determinant is 
1.18. 

The scheme given by equation (36) did not work for Test 1. Figure ^ shows the 



density and velocity obtained for Test 1; all the other schemes are able to give oscil- 
lation free solutions. There is no appreciable difference between the Frink and upwind 
reconstructions. 

Figure ([s]) shows the density for Test 2 obtained from all the three schemes. Due to the 
expansion wave, the density reaches a very low value at the middle of the computational 
domain as shown in the zoomed view in flgure (Isl-b). All the schemes are able to preserve 
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Figure 4: Solution for Test case 1 at time t = 0.2 

the positivity of density. The upwind scheme has the least overshoot of density while 
Jameson scheme has the highest. The zoomed view in figure (|5]-c) show the tip of the 
expansion wave where the solution has a corner shape. We see that the two new limited 
reconstructions proposed here give the best results in terms of preserving the corner shape. 

6.2. Transonic flow over Onera M6 wing 

This is a very standard test case for compressible flows using finite volume methods 
for which experimental results are available. The problem involves transonic flow over 
the Onera M6 wing at a free stream Mach number of 0.839 and an angle of attack of 
3.06 degrees |T]. The solution consists of shocks on the upper surface of the wing with 
a lambda structure on the wing. The tetrahedral grid used in the computations consists 
of 341797 cells. With the original vertex averaging procedure, there are 2044 vertices 
with negative weights and the minimum determinant is 1.15E-15, while with the modified 
averaging procedure, there are only 86 vertices with negative weights and the minimum 
determinant is 0.098. 

The pressure coefficient at different spanwise stations are plotted and compared with 
experimental results in figure (|6|. The Frink and upwind reconstruction schemes yield 
very similar results, while the Jameson scheme is slightly dissipative at the shocks. The 
comparison with experimental results is remarkably good for all the schemes considering 
that the grid used is not very fine. The pressure variation on the wing surface is shown 
in figure ([T]) for the three schemes with the same contour levels; the lambda shock is well 
captured by all of them. However, the Jameson reconstruction can be again seen to give 
slightly more diffused shock compared to the other two schemes. 



19 




(a) (b) (c) 

Figure 5: Solution for Test case 2 at time t — 0.15 

6.3. Supersonic flow over hemi-sphere 

This problem consists of a Mach 3 flow over a hemi-sphere which leads to the formation 
of a bow shock in front of the sphere. We use the KFVS scheme for the numerical flux 
function, which is known to be entropy consistent and gives physically correct solutions 
even for strong shocks. The radius of the hemi-sphere is one and the problem is solved 
on two grids: grid Gl has a cell size of 0.05 on the surface of the hemi-sphere with 
908330 tetrahedra while grid G2 has a cell size of 0.025 with 1934791 tetrahedra. All the 
averaging weights were positive using both methods. In the original procedure of Frink, 
the smallest determinants were 2.9e-ll and 3.5e-13. In the case of the modified averaging 
procedure, the smallest determinants for the two grids were 1.54 and 1.98 respectively. 

The variation of pressure along the stagnation stream-line is shown in figure (|8]). Along 
this line, the shock is normal to the line and we see that there are no oscillations in the 
pressure distribution. This is true of the density and velocity also but the plots are not 
shown here. The shock width is seen to decrease when we go from the coarser grid Gl 
to the finer grid G2. Figure ^ shows the pressure contours obtained on G2 using the 
three reconstruction methods. It is clear that all of them give similar results without any 
oscillations. 

6.4. Blast wave 

In order to test the limiters under strong shocks, we consider a problem with a large 
temperature/pressure discontinuity. The initial condition consists of a small region of 
high temperature. The solution develops into an expanding shock wave. In the initial 
condition, the density is uniform with a value of 1.228 while the velocity is zero. The 
temperature inside a spherical core of radius 5 meters is 8.1 x 10'' K while outside the core 
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Figure 6: Pressure coefficient for Onera M6 wing 
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Frink Jameson Upwind 



Figure 7: Pressure contours for Onera M6 wing 




(a) Grid Gl (b) Grid G2 

Figure 8: Pressure variation along stagnation stream-line for flow over hemi-sphere 
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(a) Frink (b) Upwind (c) Jameson 

Figure 9: Pressure contour for flow over hemi-sphere using grid G2 



it is 298 K. The numerical flux function is based on KFVS scheme. The computational 
domain is a cube whose sides are 81 meters and the grid consists of 3182815 tetrahedra. 
With the original averaging procedure, there was one vertex with negative weight and the 
minimum determinant was 8.7 X 10-^. With the modified averaging procedure, there was 
no negative weight and the minimum determinant was 1.5. 

The computations are performed in a time accurate manner using 3-stage Runge- 
Kutta scheme and a constant time step of 6 x 10~^ for which the CFL number is less than 
one. The Frink and Jameson reconstruction schemes lead to loss of positivity of pressure 
in the first iteration itself. Only the upwind scheme is able to give stable solutions for this 
problem. In figure ( 10 ) we compare the radial pressure variation for the first order scheme 
and the second order scheme with upwind reconstruction at two different times; note that 
the pressure jump across the shock has a large magnitude in the higher order scheme. 
It is seen that the second order scheme is also free of oscillations. Moreover, the second 
order scheme is able to resolve more variations in the solution as compared to first order 
scheme, especially at the earlier times. This indicates that the limiter is not completely 
reducing the scheme to first order. A similarity solution for the radius of the blast wave 
is given by Taylor [20] according to which the quantity R^^'^t~^ is a constant, where R 
is the radius of the blast wave. From the numerical solution of the pressure variation 
obtained with the upwind reconstructed scheme, we determine the blast wave radius and 
the results are plotted in figure (11) along with the analytical results. Except for the 
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early times, the evolution of the radius is well predicted by the current computations as 
compared to theoretical results. 

7. Conclusions 

Vertex-centroid finite volume schemes on tetrahedral grids have a simple structure 
even for second order schemes since they make use of vertex values for reconstruction. 
Maintaining the positivity of the interpolation weights is important for stability of the 
scheme. The new averaging procedure given here is shown to be more successful in main- 
taining positivity of the interpolation weights and has good properties in terms of having 
0(1) determinants. The simplified reconstruction scheme (upwind scheme) proposed here 
is shown to be theoretically more accurate as compared to the Frink scheme, and gives 
comparably similar results on test cases. By writing the scheme as a difference of vertex 
values and cell-center values, limited versions of the scheme are developed and shown to 
be stable in maximum norm for scalar conservation laws. Upwind scheme is seen to be 
less restrictive in terms of the conditions that the limiter has to satisfy. These schemes 
give oscillation free solutions when applied to Euler equations. For the blast wave prob- 
lem which has a very strong shock, only the upwind limited scheme was stable while the 
other schemes failed due to loss of positivity of density /pressure, demonstrating the good 
stability properties of the proposed reconstruction scheme. 
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Figure 11: Variation of blast radius with time. The straight hue corresponds to the relation 
i?5/2;-i =const. 



References 

[1] AGARD. Experimental database for computer program assessment. Technical Re- 
port AR-138, AGARD, 1979. 

[2] T. J. Barth. Aspects of unstructured grids and finite volume solvers for the Euler 
and Navier-Stokes equations. In Lecture Series in Computational Fluid Dynamics, 
No. 1994 -04- von Karman Institute for Fluid Dynamics, 1994. 

[3] T. J. Barth and D. C. Jespersen. The design and analysis of upwind schemes on 
unstructured meshes. Number AIAA-89-0366, 1989. 

[4] N. T. Frink. Recent progress towards a three-dimensional unstructured Navier-Stokes 
flow solver. In 32nd Aerospace Sciences Meeting, number AIAA-94-0061, Reno, NV, 
January 1994. 

[5] N. T. Frink and S. Z. Pirzadeh. Tetrahedral finite-volume solutions to the Navier- 
Stokes equations on complex configurations. Technical Report TM-1998-208961, Lan- 
gley Research Center, 1998. 

[6] C. Geuzaine and J.-F. Remade. Gmsh: A 3-d finite element mesh generator with 
built-in pre- and post-processing facilities. International Journal for Numerical Meth- 
ods m Engmeermg, 79(11):1309-1331, 2009. 



25 



[7] J. B. Goodman and R. J. LeVeque. On the accuracy of stable schemes for 2D scalar 
conservation laws. Mathematics of Computation, 45(171):15-21, 1985. 

[8] A. Harten. On a class of high resolution total-variation-stable finite-difference 
schemes. SIAM Journal on Numerical Analysis, 21(l):l-23, 1984. 

[9] D. G. Holmes and S. D. Connell. Solution of the 2D Navier-Stokes equations on un- 
structured adaptive grids. In AIAA 9th Computational Fluid Dynamics Conference, 
June 1989. 

[10] A. Jameson and J. C. Vassberg. A vcrtcx-centroid (v-c) scheme for the gas-dynamics 
equations. In N. Satofuka, editor, Computational Fluid Dynamics (2000). Springer, 
July 2000. 

[11] E. Kreyszig. Introductory Functional Analysis with Applications. Wiley, 1989. 

[12] S. H. Lui and K. Xu. Entropy analysis of kinetic flux vector splitting schemes for the 
compressible Euler equations. Zeitschrift fiir Angewandte Mathematik und Physik 
(ZAMP), 52:62-78, 2001. 10.1007/PL00001540. 

[13] J. C. Mandal and S. M. Deshpande. Kinetic flux vector splitting for Euler equations. 

Computers and Fluids, 23(2) :447 - 478, 1994. 

[14] D. J. Mavriplis. Unstructured mesh discretizations and solvers for computational 
aerodynamics. In 18th AIAA Computational Fluid Dynamics Conference, 2007. 

[15] P. L. Roe. Approximate riemann solvers, parameter vectors, and difference schemes. 
Journal of Computational Physics, 43(2) :357 - 372, 1981. 

[16] W. Rudin. Principles of Mathematical Analysis. McGraw-Hill Book Company, 1976. 

[17] D. Sharov and K. Nakahashi. Low speed preconditioning and LU-SGS scheme for 
3-D viscous flow computations. AIAA Paper 98-0614, 1998. 

[18] D. Shepard. A two-dimensional interpolation function for irregularly-spaced data. 
In Proceedings of the 1968 23rd ACM national conference, ACM '68, pages 517-524, 
New York, NY, USA, 1968. ACM. 

[19] C.-W. Shu and S. Oshcr. Efficient implementation of essentially non-oscillatory shock- 
capturing schemes. Journal of Computational Physics, 77(2) :439 - 471, 1988. 

[20] G. I. Taylor. The formation of a blast wave by a very intense explosion. I. Theoretical 
discussion. Royal Society of London Proceedings Series A, 201:159-174, Mar. 1950. 



26 



